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ABSTRACT 


By  modelling  a  picture  as  a  two-state  Markov  field,  MAP  estimation 
technqiues  are  used  to  develop  subopt iaal  but  computationally  tractable 
binary  segmentation  algorithms.  The  algorithms  are  shown  to  perform 
well  at  low  signal  to  noise  ratios,  and  analytical  procedures  are  developed 
for  estimating  the  Markov  field  transition  probabilities.  In  addition, 
extensions  of  this  approach  to  the  multi-spectral  and  multi-region  cases 
are  discussed. 


I.  INTRODUCTION 

The  image  segmentation  process  Is  a  basic  component  of  computer 
vision  systems,  and  as  such,  has  received  considerable  attention  in  the 
Image  processing  literature.  In  this  report  ve  present  a  method  for 
segmentation  of  two  Intensity  level  monochromatic  pictures  in  the  pres¬ 
ence  of  high  levels  of  additive  noise.  Although  we  concentrate  on  this 
limited  class  of  images,  our  approach  is  extendable  to  a  much  more 
general  class  of  Images  and  this  is  briefly  discussed  in  Section  VI  of 
the  paper.  In  order  to  work  with  low  signal  to  noise  ratios,  we  have 
exploited  the  spatial  dependence  of  pixel  intensity  by  modelling  the 
true  underlying  picture  as  a  Harkov  field. 

In  the  case  of  one  dimensional  digital  signals  the  Harkov  field 
model  reduces  to  a  Harkov  chain.  Harkov  chain  models  for  stochastic 
processes  have  received  considerable  attention  in  the  statistical  liter¬ 
ature  [1],  [2],  and  found  extensive  application  in  the  control,  coasnmi- 
catlon,  and  signal  processing  fields  [3]-[6].  Estimation  and  detection 
problems  associated  with  random  signals  modelled  by  Harkov  chains  can  be 
formulated  as  likelihood  maximisation  problems  where  one  wants  to  maximise 
the  joint  likelihood  of  the  data  and  the  Markov  state  sequence.  Maximisation 
of  this  joint  likelihood  is  completely  equivalent  to  generating  what  is  known 
as  the  maximum  a  posteriori  probability  estimate  of  the  Markov  chain  [18] • 

In  the  one  dimensional  case,  this  leads  to  etegant  and  reasonably  efficient 
dynamic  programing  algorithms  for  computing  the  estimate  which  maximises 
the  joint  maximum  a  posteriori  probability  of  the  entire  data  string  (digital 
signal)  in  a  sequential  manner.  Unfortunately,  this  approach  does  not  gener¬ 
alise  in  a  natural  way  to  the  case  of  two  dimensional  signals  such  as  digital 
images.  As  a  result,  this  has  somewhat  hampered  the  use  of  MAP  formulations 
in  image  processing.  The  few  exceptions  [7]-[12]  are  discussed  more  carefully 


below. 
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Recent  work  by  Kaufman  and  co-workers ,  [7],  and  Therrien  [8],  [9] 
have  also  made  use  of  Markov  field  models  and  MAP  formulations.  In  [7]  this 
approach  was  combined  with  reduced  update  Kalman  filtering  techniques 
for  image  enhancement,  while  in  [8]  and  [9]  it  was  combined  with  two- 
dimensional  autoregressive  texture  models  for  texture  based  segmentation. 
However,  the  algorithm  described  in  [7]  and  one  of  the  algorltham 
described  in  [8]  and  (9]  fall  to  exploit  the  true  spatial  dependence 
Imposed  by  the  model.  Zn  particular  they  don't  attempt  to  maximise  a 
joint  likelihood  of  all  the  data  but  rather  the  individual  likelihoods 
at  points  or  In  small  regions  within  the  image.  This  problem  is  partially 
overcome  by  a  second  multi-pass  or  Iterative  algorithm  proposed  in  [8] 
and  [9].  However,  the  relation  between  the  estimate  obtained  using  this 
approach  and  a  true  joint  MAP  estimate  is  unclear. 

Cooper  and  Elliott  [10]-[12]  have  applied  MAP  techniques  to 

boundary  estimation  in  noisy  Images  as  well.  Although  boundary  estimation 
can  be  formulated  as  a  one-dimensional  signal  estimation  problem  where 

the  independent  variable  la  arc  length  along  the  boundary,  it  is  Inter¬ 
esting  to  note  that  even  in  this  case  the  two  dimensional  Image  data 

necessitated  the  use  of  a  suboptlmal  algorithm. 

In  this  report,  we  have  taken  a  very  simple  Markov  field  model,  and 

have  developed  an  algorithm  which  approximates  the  behavior  of  an  optimal 
sequential  estimation  algorithm.  It  makes  use  of  two  stages  of  dynamic 
programming.  In  the  first  stage  a  "generalized"  dynamic  programming 
algorithm  is  applied  to  each  row  of  the  image.  This  yields  a  set  of 
candidate  segmentations  for  each  row.  In  the  second  stage,  a  final  seg¬ 
mentation  is  "pieced"  together  from  the  candidate  row  segmentations  using 
dynamic  programming  as  well.  The  algorithm  requires  only  a  single  pass 
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over  the  image  data,  and  the  versioh  for  moderate  signal  to  noise  ratios 
can  be  performed  in  a  highly  parallel  fashion.  The  version  which  works 
best  at  low  signal  to  noise  ratios  requires  a  sequential  raster  proces¬ 
sing  of  the  luge. 

The  paper  is  organized  as  follows.  In  Section  II  the  problem  of 
interest  is  formulated  and  the  image  model  is  discussed.  In  Section  III 
the  basic  suboptimal  algorithm  is  derived,  while  in  Section  IV,  analysis 
tools  are  developed  to  allow  estimation  of  the  Markov  field  transition 
probabilities.  Section  V  presents  some  modifications  to  the  original 
algorithm  to  improve  its  performance  at  low  signal  to  noise  ratios,  and 
extensions  Of  this  approach  to  a  larger  more  Interesting  class  of  Images 
is  given  in  Section  VI.  Finally,  some  concluding  remarks  are  made  in 
Section  VII. 


IX.  MODELING  ASSUMPTIONS  AND  PROBLEM  FORMULATION 


Let  a  digitized  monochromatic  picture  be  modeled  as  an  (N^zN^)  matrix 


B  with  components,  b^,  representing  the  Intensity  or  grey  level  of  pixel 


(l,j).  Although  this  model  will  be  generalized  In  Section  VI ,  Initially 
b^  will  be  restricted  to  take  on  only  one  of  two  values,  r^  and  ,  r^  > 
.  Thus  the  pixels  in  an  Image  can  be  classified  into  two  sets  and 


S2  where 


SL  *  Ui,j):  bij"riJ 
S2-{(i,j):  b±j-r2} 


(la) 

(lb) 


It  will  be  assumed  that  an  observed  Image,  G,  represents  a  picture  B 
corrupted  by  a  stationary  random  noise  field  W  so  that  G-B+W  or  equiva¬ 


lently,  The  random  variables  w^  are  assumed  to  be  lndepen- 


2  2 
dent,  Gaussian,  zero  mean,  and  of  known  variance  0  ,  denoted  wjj ~N(0,o  ). 


In  this  report  we  consider  the  problem  of  segmenting  the  Image,  or 
equivalently,  generating  estimates  and  §2  of  the  sets  and  S2<  In 
the  context  of  the  model  for  B,  this  problem  can  be  phrased  as  a  two  set 
classification  problem  where  pixel  (i, j)  is  assigned  to  or  §2  using 


the  noisy  observation  g^.  Since  the  conditional  probability  distribu¬ 


tions 


P(g1j<t|(i,J)eS^-  /  (2wo2)_1^2exp(-(x-rk)2/2o2)dx 


(2) 


k*l,2 

are  assumed  known,  maximum  likelihood  classification  techniques  could  be 
employed.  This  would  result  in  a  simple  thresholding  algorithm  with 


pixel  (i,j)  assigned  to  Sj  if  g^A  and  to  §2  if  where  A^rj+tj)^ 

This  method  only  works  well  if  the  signal  to  noise  ratio  S/N  A  A la  >  2, 


where  A  A  rj-r2> 


****** 


i 
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s. 


In  order  to  derive  e  more  robust  segmentat Ion  procedure,  it  is 
necessary  to  inpose  additional  structure  on  B.  Pixels  in  the  set 
will  be  sssuaed  to  appear  in  connected  subsets,  or  blobs,  as  Illustrated 
in  Figure  1.  This  structure  can  be  modeled  by  a  Markov  field.  Assuming 
the  support  for  this  field  to  be  limited  to  a  pixels  four  nearest  neigh¬ 
bors,  this  Markov  process  can  be  characterized  by  the  transition  proba¬ 
bilities: 


P(bij“rklB)  "  P(bij"rklbi-l,J*bi+l,J‘bi,J-l»bi,j+l) 
&  Pijk  * 


(3) 


This  model  is  fairly  general  in  nature,  is  direction  invariant,  and  can 
be  used  to  model  a  large  class  of  clusters.  Unfortunately,  it  is  diffi¬ 
cult  to  identify  appropriate  transition  probabilities.  In  addition, 
optimal  segmentation  based  upon  such  a  model  would  be  computationally 
impractical.  Thus  in  the  section  to  follow,  we  introduce  a  number  of 
restrictions  and  approximations  in  order  to  arrive  at  a  suboptlmal  but 
computational  tractable  segmentation  algorithm  based  upon  dynamic  pro¬ 
gramming  concepts. 
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III.  SUBOPTIMAL  SEGMENTATION 


In  order  to  motivate  our  subopt  imal  algorithm  ve  first  formulate  an 
optimal  algorithm  which  maximises  the  joint  likelihood  of  the  Image  data, 
and  the  state  of  the  Markov  field  assuming  the  transition  probability 
structure  given  by  (3). 

An  Optimal  MAE  Formulation 

Let  S  ”  represent  an  estimate  of  S  ■  {S^.S^},  and  let  *.(•) 

represent  a  log-likelihood  function.  The  joint  log-likelihood  of  S  and 
the  observed  data  G  can  be  written  as: 

i(§,G)  »  t(G|§)  +  i(S)  .  (A) 

A 

The  log-likelihood  of  the  data  conditioned  on  the  segmentation  S  is 


£<G|S)  -  C  -  I  .  “V  "  1 

(i,j)eS1  2a  13  1  (i,j)eS2  2o 


■h  (*ii-r2)2  (5) 


where  C  is  a  constant  Independent  of  the  estimated  segmentation  S.  The 
log-likelihood  of  the  estimated  segmentation  is 


\ 


£(S)  -  I  A  to  P  +  1 

(i,j)eS1  U.j)eS2 

Combining  (A) ,  (5) ,  and  (6)  yields 


to  P 


ij2 


(6) 


i(S,G)  »  Z  (toP  -  “  +  1  »  (toPH2  ‘  “T  (8ii“r2) 

(i.J)^  131  2o2  i3  1  (i,j)eS2  1JZ  2oZ  13  2 

(7) 


For  optimal  segmentation  (7)  would  have  to  be  maxi¬ 

mised  over  all  possible  segmentations  S.  Even  if  appropriate  transition 
probabilities  could  be  identified,  maximisation  of  (7)  is  computationally 
impractical  since  it  cannot  be  maximised  in  a  sequential  manner.  In 
order  to  derive  a  computationally  feasible  algorithm  i(G,§)  will  be 
approximated  by  a  function  l(G,§)  which  can  be  maximised  by  processing 
the  observation  matrix  G  in  a  sequential  manner.  It  might  be  noted  that 
the  log  likelihood  (A)  differs  from  the  likelihood  traditionally  used  for 
MAP  estimation  by  only  a  term  which  is  Independent  of  §.  For  the  lattei;  one 
would  use  l(a/G)  •  i(G/C)  +  I(S)  -  1(G). 
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A  Subopt  Inal  Likelihood  Function 

To  derive  the  approximation  t(G,S)  we  will  restrict  the  Kerkovisn 
dependence  of  to  one  side.  In  particular  we  will  assume 

Pijk  "  P(bij"rklbi-l,J,bi,j-l)  * 

Use  of  (8)  in  (7)  would  still  require  sequencing  through  the  image  in 

two  directions  simultaneously  for  maximisation.  Since  no  efficient 

algorithm  has  been  formulated  for  such  a  maximisation  an  additional 

simplifying  approximation  will  be  made.  In  particular,  we  will  approxi¬ 


mate  PiJk  by  P  k  where 

Pijk  "  Rijk  Cijk 

Cijk  "  P(biJ  "  rklbi-l,J* 

V  “  P(biJ  “  rklbi,j-l> 


(9a) 

(9b) 

(9c) 


Ry k  represents  a  transition  probability  from  the  preceding  pixel  in 
the  same  row  and  C^k  represents  a  transition  probability  from  the  pre¬ 
ceding  pixel  in  the  same  column.  In  the  next  section  we  will  present 
some  methods  for  estimating  R^k  and  C^k  from  a  priori,  knowledge  of 
Image  properties.  It  might  be  noted,  however,  that  if  one  starts  with 
knowledge  of  the  probabilities  P^k>  R^k  and  C^k  can  be  calculated 
recursively.  Using  (9)  we  can  define  £(G,S)  by 


£(G,S) 

*  Hr(g,s)  + 

(10a) 

iR(G,S) 

-  1 

(i.DeSj 

(inR^ 

1 

2o2 

(10b) 

+  £  . 
(i,j)eS2 

<W*lj2 

1 

_2o2 

(«lTr2)2) 

jL(S) 

-  Z 

inC. 

+  Z 

taC _ 

(10c) 

C 

(i.DeSi 

1J1 

(1,3)  eS2  *** 

(10c) 
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Thus  by  using  (9)  ws  have  decomposed  our  likelihood  into  two  components. 
The  first,  iR,  contains  the  Image  data  and  the  tow  transition  probabili- 

the  second,  ,  contains  only  the  coluan  transition  probabili¬ 
ties. 

Subopt Inal  Algorithm  for  Maximisation  of  l(G.S 


In  order  to  use  computationally  efficient  sequential  algorithms  for 
estimating  S  it  will  be  necessary  to  maximize  the  approximate  likelihood 
(10)  in  a  subopt laal  manner.  In  order  to  explain  the  algorltlm  It  will 
be  helpful  to  introduce  some  additional  notation.  Let 
G(i,J)  A  first  J  elements  in  lth  row  of  G 
S(i, J)  -  {E1(l,J),82(i,J)) 

"  subset  of  S  containing  first  J  pixels  in  row  i  of  image. 

We  can  then  write 


S  -  U  S<i,lO 
1-1  1 


I  -M 


^(G'S)  -  I  ^(Ga.NjMU,^))  (11a) 

l  (G(i,N  ),S(i,S2))  -  Z  .  (jtaR  1  (g  r  )2} 

*  2  (1,  j)  eS^i,^)  1J1  2o2  y  1 

+  1  (taR  -  —^2  (g  -r  )2)  (lib) 

(i,j)ES2(i,M2)  1JZ  2oZ  iJ  2 

Our  algor lthn  consists  of  two  stages. 

Stage  1 ;  Por  each  row  1,  1  <  1  <  Hj,  we  calculate  the  M  row  sets 

S*(i»»2^»  1  1  M, which  yield  the  M  largest  values  of  fR(C(i,»2),S(i,»2)), 

Observe  from  (11)  that 

iR(C(i,J),S(i,J))  -  1 (G(i,J-l),S(i,J-l))  +  £(g, r)  (12) 
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a  12 

where  *(gjj)  ”  - 2  ^8u“rk^  *  *or  or  depending  upon 

2a 

pixel  (1,J)  is  assigned  to  set  or  S2«  Recursion  (12)  In  con¬ 
junction  with  the  principle  of  optlaality  [13],  Imply  the  existence  of  a 
forward  dynamic  programming  algorithm  for  finding  the  estimate  <W> 
which  maximizes  (lib)  •  This  type  of  dynamic  programming  algorithm  is  extended 
in  a  straight  forward  manner  to  generate  the  M  most  likely  sets  Sm(i,N2), 

1  £  m  _<  M.  A  detailed  description  of  the  algorithm  can  be  found  in  Appendix  A. 

Stage  2:  Given  the  sets  Sm(i,N2>,  1  £  If  1  <  1  <  M,  we  calculate 
our  complete  estimate  S  as 


Ni 

1  -wi 

§  -  U  S  1(i,N  ) 

1-1  * 

N 

where  the  sequence  of  sets  {S  (i.N^) maximize 


N 


1  - 


£<G.S)  -  Ji[lR(G(i,M2),S“(i,N2))  +  ic(sm(i,N2))] 


(13) 


Define 


£(G,S,I) 


I 

S  iR(G(ifN2),s”(i,N2))  +  ic(S*(i,N2)) 


(14) 


so  that  £(G,S)  “  £(G,S,Nj).  By  observing  that 

£(G,§,1)  -  £(G,§,I-1)  +  IR(G(I,N2),Sm(I,N2)> 

+  £c(§“(I,N2))  (15) 

He  can  again  apply  the  principle  of  optimality  to  develop  a  forward 

N1 

dynamic  programming  algorithm  for  generating  {S  (i,N.) }  .  The  details 

1  i-1 

of  this  algorithm  can  also  be  found  in  Appendix  A.  At  this  point  we  make 
the  following  remarks  with  respect  to  this  algorithm. 

Remark  1:  In  developing  this  algorithm  we  have  traded  off  optimality  in 
order  to  reduce  one  two-dimensional  maximization  problem  into  (Nj+1)  one- 
dimensional  problems.  To  this  point  in  our  experimentation  we  have  chosen 
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I 

*2 

M  such  that  1  £  M  <_  8.  Actually  as  M  approaches  2  the  algorithm  converges 
to  an  optimal  algorithm  in  the  sense  that  it  generates  an  estimate  S  which 
maximizes  (10). 

Remark  2 ;  This  algorithm  is  attractive  computationally  since  the  time 
consuming  Stage  1  is  ideally  suited  £or  parallel  Implementation,  where 
processors  simultaneously  process  each  row.  Furthermore  the  amount  of 
memory  needed  to  store  the  Stage  1  results  needed  for  execution  of  Stage  2 
are  modest.  Each  partition  Sm(i,N2)  can  be  stored  as  a  simple  bit  string 
of  length  where  a  one  in  position  j  Implies  pixel  (l,j)  is  assigned  to 

a  A 

set  and  a  zero  implies  pixel  (i,j)  is  assigned  to  set  S^.  Thus  for 
each  row  of  the  image  we  need  to  store  M  bit  strings  of  length  and  M 
numbers  corresponding  to  the  value  of  £R(G(i,N2),  Sm(i,N2)). 

Figures  2  and  3  give  examples  of  the  algorithms  performance  in  finding 
an  ellipse  imbedded  in  noise  such  that  (r^-r2>/o  ^  A/o  *  2  and  1  respec¬ 
tively.  The  results  are  compared  with  a  simple  thresholding  algorithm. 

For  these  simulations  we  chose  !^1,4,  and  8,  and  we  used  the  transition 
probabilities 


T  if  pixel  <i,j-l)e§k 
1-T  otherwise 


(16a) 


T  if  pixel  (i-l,j)eSk 

1-T  otherwise 


(16b) 


where  T».95.  This  choice  for  R^k  and  C^k  imply  both  stationarity  and 
symmetry  in  the  1  and  j  directions.  Methods  for  choosing  T  are  discussed 
in  the  following  section.  Observe  that  the  algorithm  as  given  above,  is 
subject  to  burst  type  errors  where  strings  of  pixels  are  mlsclassified. 
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This  is  •  coboo  ittrlbota  of  Algorithms  which  estimate  Markovian  pro¬ 
cesses  using  ■erlsBS  likelihood  techniques.  Fortunately  these  type  errors 
are  easily  coepen sated  for  using  simple  post-processing  algorithms.  This 
and  other  Improvements  to  the  basic  algorithm  are  discussed  in  detail  in 
Section  V. 
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IV.  ALGORITHM  ANALYSIS  FOR  PARAMETER  CALCULATIONS 
In  this  section  a  simplified  error  analysis  is  presented.  Using 
techniques  similar  to  those  used  in  [4],  we  calculate  the  probability 
that  a  given  Incorrect  row  segmentation  is  more  likely  than  the  under¬ 
lying  true  row  segmentation.  This  probability  is  then  used  to  obtain  an 
estimate  of  the  parameter  T  in  terms  of  certain  a  priori  information 
about  the  structure  of  B,  Some  examples  of  the  algor itbms  performance 
with  various  values  of  T  will  then  be  presented. 


Simplified  Error  Analysis 

We  "»»<»»  the  same  assumptions  about  the  transition  probabilities 

and  R  .  that  were  used  for  the  simulations  described  in  Section  III.  In 
ijk 

particular,  we  assume  statlonarity  and  symmetry  so  that  (16)  holds.  Let 
S(i,N2)  -  {SjU.N^.SjU.Nj)}  be  an  Incorrect  estimate  of  the  true  seg¬ 
mentation  S(i,N2)  -  (S1(i,N2)>S2(i,N2)}  of  row  i.  Consider  the  probability 
PE(i),  that  S(i,N2)  is  more  likely  than  S(i,N2>.  Using  the  stage  1  like¬ 
lihood  function,  *(•.•) »  we  obtain 

PE(i)  -  P(£RG(i,N2) ,S(i,N2))  >  £R(G(l,N2),S(i,N2)))  (17] 


where 


P(Q(i)  >  0) 


Q(i)  ^  ER(G(i,N2),S(i,N2))  -  £R(G(i,N2),S(i,N2)) 

Using  (lib)  and  expanding  the  quadratic  terms,  Q(l)  can  be  rewritten 


Q(i)  - 


l  (in  R  +  (2g  r.-r^)) 

(i,j)eS1(i,N2)  131  2<j  13  1  1 


+  £  (in  R.,,  +  -K  (2g,.r,-r?)) 

(i,j)e$2(i,N2)  132  2o2  13  2  2 


(i,j)eS1(i,N2) 


(in  Rijl  +  2  <28ijri“r1>3 


In  these  expressions  R^k  will  take  the  value  T  if  pixels  (i,j-l) 
and  (i,j)  are  placed  into  the  same  set,  and  will  take  the  value  (1-T) 
otherwise.  Let  t(i)  be  the  number  of  times  pixels  (i,j-l)  and  (i,j) 
are  classified  differently  in  S (i,^) ,  and  let  t(i)  be  defined  analogously 
for  S(i,N2).  If  §^(i,N2)  is  defined  as  the  set  of  pixels  contained  in 
Sk(i,N2)  and  not  in  Sk(i,N2),  i.e.,  §k(i,N2)  -  Sk(i,N2)  -  Sk(i,N2)  then 
Q(i)  can  be  simplified  to 


Q(i) 


L  (^-^^(r.-r-Jg..)  + 

(i,j)e§2(i,N2)  12  3  (19) 

,,  | _ .  (t2-rr2(r2-ri)8i))1+<t<1)  -  ‘<im" 

(it J) eS^(i,N2) 
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for  £(1,^).  Let  n^(l)  be  the  number  of  pixels  in  the  set  S^d,^)* 

Then,  conditioned  on  the  knowledge  thst  Sd.N^)  is  correct 

a2  a  ? 

Q(i)  -  (51(i)  +  52(i))  +  (t(i)  -t(i))tn(T/l-T),  (f)  («jd)  +  52CD) 

(20) 

Thus 


pE<i) 


0D  2 

/  (2ir)”1  2exp(-  ^-)ds 

o(i)  2 


(2le) 


where 


a(i)  - 


(51(i)-hl  (i))1/2 


_  ,1,  _ 


(n^(i)+n2(i)) 


1/2  in^l-T^ 


(21b) 


As  expected,  Pg(i)  goes  to  zero  as  the  signal  to  noise  ratio, 

A/o,  goes  to  infinity.  The  contribution  of  the  transition  probabilities 
to  Pg(i)  diminishes  with  Increasing  signal  to  noise  ratio,  and  the  tran¬ 
sition  probability  term  only  influences  Pg(i)  when  t(i)l*t(i) ,  or  when 
the  number  of  state  changes  in  S(i,N2)  is  different  than  the  number  in 
S(i,N2). 

Harkov  Parameter  Determination 

Using  (20)  and  (21),  bounds  on  the  parameter  T  can  be  deter¬ 
mined.  These  bounds  can  be  arrived  at  by  looking  at  two 
limiting  cases.  One  case  is  when  T  is  very  large  so  that 
the  probability  of  changing  states  is  low.  Here  it  would  be  expected 
that  an  object,  a  connected  set  of  pixels  in  either  or  S2,  might  be 
deleted  from  the  estimated  segmentation,  §,  due  to  the  low  probabil¬ 

ity  of  any  state  translations  from  to  S2  or  S2  to  S^.  The  other 
extreme  is  T  very  smell  which  makes  the  probability  of  changing 


states  larger.  In  this  case  s  number  of  extraneous  transitions  might 


14 


be  Included  in  5.  This  situation  produces  what  might  be  called 'jfalse 
object^',  and  will  be  considered  first. 

For  the  first  case,  let  3(1, M^)  be  such  that  all  pixels  In  a  connected 
set,  rx,  in  row  i  are  misclassified.  Furthermore,  it  will  be  assumed 
that  the  set  is  such  that  the  mlsclassif icatlon  causes  Sd.N^)  to  have 
2  state  transitions  more  than  Sd.M^).  An  illustration  ef  such  an  occurrence 
appears  in  Figure  4a. 

If  the  set  I^ls  n^  pixels  long,  then  the  probability  that  3(1, N^) 
will  be  more  likely  than  S(i,N2>  is  given  by  (17)  with 
(n  )1/2 

a(i)  -  ~  l2  (|)  +  (|)  ~~^J2  4n(T/l-T)  (22) 

“l 

Thus  to  assure  that  P  (i)  is  small,  l.e.  P_(i)  <  e,  T  should  be  picked 

£  £ 

such  that  a(i)  is  greater  than  or  equal  to  some  value  6(c).  From  (22  ) 


this  is  equivalent  to  requiring 

T  >  exp(8(e)x-x2)  (l+exp(B(e)x-x2)) 


(23a) 


where 


.H 

|A, 

2  >o 


(23b) 


Mote  that  the  left  hand  side  of  this  inequality  is  maximised  at 
x  »  p(e)/2,  thus 

T  >  exp  /  ( 1+exp  (~^^-) ) 

will  guarantee  that 


P_(i)  <  /  (2ir)"1/2  exp(-s2/2)ds  -  e  (25) 

E  8(e) 

For  example,  if  we  set  e*.01  then  8(e)  *  2.33.  This  will  require  that 
T  >_  .80  which  is  quite  reasonable. 


If  (22)  is  minimised  with  respect  to  n^,  we  see  that  given 


(~)and  T,  the  most  probable  length  of  a  false  object  la  approximately 
n*  -  4fn(T/l-T)(7-)2  (26] 

A 

Also,  it  should  be  noted  that  the  event  of  the  row  segmentation 


S(1,N2) 


being  more  likely  than  S(1,N2>  la 


dependent  only  on  the  pixels  in  1^.  Therefore  one  estimate  for  the 

nuaber  of  possible  Independent  sites  for  such  a  set  to  occur  might  be 

(Hi  N2)/(nl*fl) Since  the  corrupting  noise  was  assumed  white,  the 

number  of  such  sets,  in  the  segmentation,  S,  can  be  similarly  estimated  as 

ni  ”  niVe  ^*+1)  (27) 

In  practice,  however,  this  estimate  seems  to  be  consistently  low  by 
a  factor  of  approximately  four. 

Thus,  given  an  upper  bound  on  the  probability  of  S  containing  a 
false  object  in  any  given  position,  a  lower  bound  on  T  can  be  obtained. 

Once  T  is  chosen  both  the  most  probable  length  of,  and  the  expected 
number  of  such  sets  can  be  estimated  by  n*  and  4p^  respectively. 

In  order  to  obtain  an  upper  bound  on  T  the  second  case  must  be 
considered.  Here,  it  is  assumed  that  an  entire  connected  set  of  pixels, 
r2.  In  S^(1,N2),  k^l,2,  has  been  falsely  classified  in  S^d.Nj).  By 

entire,  we  mean  that  the  pixels  immediately  to  either  side  of  r  2  are 
not  elements  of  S^(i,N2).  In  thiB  case  §(1,N2)  has  two  fewer  state 
transitions  than  S(1,N2>.  This  case  is  illustrated  in  Figure  4b. 

Again,  it  is  assumed  that  §(i,N2)  includes  no  other  errors  and 
that  is  n2  pixels  long.  In  this  case  the  probability  of  an  error  Pg(i), 

+ 

n*+l  is  used  rather  than  n*  to  insure  each  connected  set  is  separated 
from  the  others  by  at  least  one  pixel.  Thus,  each  "false  object" 
will  contribute  2  extraneous  translations  to  t(l). 


VUi 


•s  defined  by  (17)  end  (18) ,  is  given  by  the  Integral  (21a)  where 


a(i) 


(28) 


If  we  use  similar  arguments  to  those  of  the  previous  case.  In  order  for 
Pg(i)  e>  we  would  require  a(l)  >_  8(e).  Using  (28)  in  this  case  yield 
the  upper  bound  on  T: 

jj  jj 

T  <  exp(^  (j)2  -  8(c)  ^§-  (A))/l+exp(^(A)2  _  B(e)  (A))  (29) 

Minimization  of  (29)  with  respect  to  n.  and  (A)  yields 

z  o 

2  2 

T  <  exp(-  -^-)/(l+exp(-  -^-))  .  (30) 

This  Implies  T<. 5  which  is  inconsistent  with  the  constraints  Imposed 
on  T  by  the  lower  bound  (24)  derived  from  analysis  of  case  1.  More 
Importantly,  it  Is  Inconsistent  with  our  basic  picture  assumption  of 
relatively  large  region  clusters. 

This  problem  can  be  overcome  by  assuming  specific  values  for  the 
signal  to  noise  ratio,  (A/a) ,  and  for  the  length  n^  of  the  misclasslf led 
pixel  string  r^.  In  this  case  if  T  is  chosen  to  satisfy  (29)  for  a 
specific  e,  then  any  string  of  length  greater  than  n^  would  have  lower 
error  probability.  Tables  1  and  2  give  the  upperbound  on  T  as  defined 
by  (29)  for  values  of  e  ranging  from  .16  to  .02  and  values  on  n^  ranging 
from  4  to  40.  In  Table  1  A/o»l  and  in  Table  2  A/o*2. 

Figure  5a  shows  a  test  Image  containing  ellipses,  triangles  and 
circles  with  different  orientations  and  sizes,  and  Figure  5b  shows  this 
same  Image  with  additive  Gaussian  noise  such  that  A/o  •  1.  The  results 


of  applying  our  algorithm  with  T-.5,  .8,  .9,  and  .95  are  shown  in 
Figures  5c-5f.  Mote  that  for  small  values  of  T  many  additional  state 


transitions  occur  and  false  objects  appear.  On  the  other  hand,  large 
values  of  T  tend  to  suppress  state  transitions  causing  case  2  situations. 
Figures  6a-6c  show  the  algorithm  performance  at  A /a  m  2  and  T  -  .95  and 
.5.  With  T  »  .5  our  algorithm  reduces  to  maximum  likelihood  thresholding 
as  discussed  in  Section  II. 
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Table  1 .  Upperbounds  on  T  determined  from  Equation  (29) 
with  A/o  «  1. 


Table  2.  Upperbounds  on  T  determined  from  Equation  (29) 
with  A/o  ■  2. 


V.  MODIFICATIONS 


As  is  apparent  from  our  previous  examples,  the  basic  algorithm  as 
described  in  Section  III  performs  well  for  signal  to  noise  ratios  of  1.5 
or  greater.  However,  its  performance  is  disappointing  at  signal  to 
noise  ratios  near  or  below  1.0.  An  obvious  approach  to  Improving  per¬ 
formance  would  be  to  Increase  M,  the  number  of  row  alternatives  stored 
during  stage  1,  to  a  value  greater  than  Unfortunately,  as  is  dis¬ 

cussed  below,  this  is  impractical.  Consequently,  we  present  two  alterna¬ 
tive,  but  more  practical  modifications,  which  when  applied  jointly  yield 
reasonable  performance  at  signal  to  noise  ratios  of  1  or  below. 

Effects  of  Increasing  M 

As  can  be  seen  In  Figures  2  and  3,  as  M  Increases  algorithm  perfor¬ 
mance  does  improve.  However,  certain  row  estimates  experience  little  or 
no  Improvement.  This  is  especially  true  in  the  case  A/a  "  1.  To  under¬ 
stand  the  effect  of  Increasing  M  one  must  realize  that  stage  2  of  the 
algorithm  can  not  generate  a  reasonable  picture  segmentation  if  the  true 
row  segmentation,  or  a  reasonable  estimate,  is  not  contained  In  the  set 
of  candidate  row  cstlsmtes  generated  during  stage  1.  With  this  in  mind 
we  considered  three  troublesome  rows  In  Figure  3c,  namely  rows  11,  27, 
and  32,  extracted  them  from  the  image,  and  applied  the  stage  1  algorithm 
to  them  with  M"64.  The  20th  and  12th  most  likely  alternatives  for  rows 
11  and  13  respectively  did  yield  an  Improved  estimate  of  their  true 
segmentations.  However,  no  such  alternative  appeared  for  row  27.  This 
Implies  that  Improved  performance  is  only  possible  if  M  is  increased 
substantially.  This  in  turn  would  lead  to  Impractical  time  and  storage 
requirements.  As  a  result  we  have  considered  two  other  alternative 
modifications.  The  first  is  a  simple  post-filtering  operation.  The 


second  Involves  a  modification  of  stage  1  of  the  algor it  tan. 
Post-Filtering 


The  errors  incurred  by  the  algorithm  are  generally  "burst"  type 
errors  where  a  sequence  of  pixels  in  a  given  row  are  mlsclassified.  Thus, 
one  tray  to  remove  the  effects  of  such  errors  is  to  post -filter  the  re¬ 
sulting  binary  Image  representing  the  segmentation.  There  are  many  types 
of  filters  which  could  be  used,  and  we  show  the  results  of  applying  one 
simple  type.  Ve  reassign  the  state  of  a  pixel  if  the  majority  of  the 
pixels  in  a  5x1  window  centered  on  the  pixel  are  in  a  different  state. 

This  is  just  a  special  case  of  a  median  filter.  The  result  of  applying 
this  filter  to  the  binary  segmentation  shorn  in  Figun  5f  is  given  in 
Figure  7a,  while  Figure  7b  shows  the  results  of  applying  the  filter  to 
Figure  6c.  For  purposes  of  comparison  we  have  also  filtered  the  segmentations 
shown  in  5c  and  6b  which  were  obtained  by  simple  thresholding.  Since  the 
errors  in  these  images  do  not  exhibit  any  directional  bias,  we  used  a  (3x3) 
window.  As  can  be  seen  by  comparing  Figures  7c  and  7d  to  5c  and  6b,  respec¬ 
tively,  the  approximate  MAP  algorithm  produces  significant  Improvement  at 
A/cr-  1  and  only  minor  improvement  at  A/o  »  2.  Thus  the  additional  complexity 
of  the  MAP  algorithm  can  only  be  justified  when  working  with  signal  to  noise 
ratios  lower  than  2. 

Non-Parallel  Algorithm 

Although  post-filtering  produces  significant  improvement  in  perfor¬ 
mance  when  compared  with  the  output  of  the  basic  algorithm,  one  weald 
still  prefer  improving  the  performance  of  the  basic  algorithm  before 
post -processing.  In  view  of  the  discussion  above,  our  goal  is  to  improve 
the  list  of  candidate  row  estimates  generated  by  stage  1.  This  can  be 


done  by  using  information  from  the  previously  processed  rove.  Unfortun¬ 
ately,  to  do  this  we  must  remove  the  parallelism  from  the  algorithm  and 
process  the  rows  sequentially. 

In  this  modified  approach  it  makes  more  sense  to  interweave  stages  1 
and  2  of  the  algorithm.  That  is,  after  processing  row  1  with  a  modified 
version  of  the  stage  1  algorithm  to  obtain  the  most  likely  candidate 
segmentations,  we  then  perform  one  step  in  the  stage  2  dynamic  programming 
algorithm.  At  this  point  we  have  M  distinct  binary  segmentations  stored 
for  the  sub- image  consisting  of  the  first  1  rows. 
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The  first  stage  algor lthm  is  modified  to  make  use  of  the  dynamic 
programming  solution  for  the  previous  i-1  rows  as  well.  In  particular 
we  now  find  the  M  candidate  estimates  which  yield  the  largest  values  of 
the  modified  row  likelihood 


iMR(G(i,M2),S(i,N2))  -  lR(G(i,N2)S(i,N2)) 


+  E  in  H  + 

(i,j)e§1(i,U2)  13 


E 

(i,j)e§2(i,N2) 


to  H 


iJ2 


(31) 


where  iR(0  is  given  by  (lib),  and 


HU- 


R 


if  biJ  "  bi-l,j 


1-Tr  otherwise 


(32) 


A  A 

The  term  b  is  the  state  estimate  for  pixel  (l,j)  and  b.  1  is  the 
most  likely  state  estimate  for  pixel  (i-1, j)  after  running  both  stages 
of  the  algorithm  on  the  first  i-1  rows  of  the  image.  In  this  case  we 
have  Incorporated  information  about  our  segmentation  of  the  previous 
rows  in  making  decisions  regarding  new  rows.  In  particular  we  penalize 
state  changes  in  the  vertical  direction  as  well  as  the  horizontal  direc¬ 
tion  when  performing  our  initial  row  segmentations.  Observe  that  the 
addition  of  the  terms  in  (31)  does  not  change  its  recursive  nature 
and  hence  dynamic  programming  can  still  be  used  to  maximize  (31). 

Finally  we  point  out  that  before  performing  the  stage  2  algorithm 
on  a  given  row  the  terms  are  subtracted  out  of  the  row  likelihoods, 
£mr(')>  reducing  it  back  to  the  form  of  &R(’)  given  in  (lib).  This 
assures  that  the  final  Image  likelihood  is  consistent  with  our  Markov 
field  model. 

Estimating  T  and  TR 


In  order  to  Implement  this  modified  algorithm  we  must  estimate 
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P*  *  y1PE1  +  y2PE2  +  y3PE3 


appropriate  values  for  T  and  TR.  The  approach  taken  is  similar  to  that 

used  for  the  original  algorithm  and  the  details  can  be  found  in  [11]. 

However,  we  briefly  outline  the  necessary  changes. 

To  estimate  T  alone  we  considered  two  types  of  row  errors.  To 

simultaneously  estimate  the  two  parameters  T  and  T  we  consider  three 

K 

error  types.  Case  1  is  identical  to  the  first  "false  object"  case  used 
in  the  earlier  analysis.  Case  2  corresponds  to  the  probability  of 
missing  the  top  row  of  an  object  while  Case  3  consists  of  the  probability 
of  a  protusion  of  one  pixel  width  extending  from  an  object.  If  P^,  P£2 » 
and  Pg^  are  defined  as  the  error  probabilities  for  each  of  these  cases 
reasonable  values  for  T  and  TR  can  be  determined  by  minimizing 

(33) 

The  positive  parameters  Y^»  Y2»  and  are  adjusted  depending  upon  what 
type  of  error  is  more  reasonable  in  a  given  application.  In  our  case 
we  were  most  concerned  with  eliminating  "false  objects." 

Thus  we  used  Yj“100  and  Y2“Y-j  *  1*  For  A/o  *  0.8,  minimization  of  (33) 
yielded  T-.83,  TR  -  .55,  P£1  -  6.6xl0_4,  P£2  -  .81  and  P£3  -  .07. 

Interestingly,  as  one  would  expect  with  this  algorithm,  PR2  is  large. 

That  is  there  is  a  high  probability  of  not  estimating  the  first  row  of  objects. 
Figure  8a  shows  the  image  of  Figure  5a  with  additive  noise  such  that 
A/o  *  .8  while  Figure  8b  shows  the  result  of  applying  the  modified  algo¬ 
rithm  to  this  image  with  T  and  TR  as  determined  above.  Figure  8c  is  the 
post  filtered  version  of  Figure  8b.  Figure  9a  shows  the  result  of  apply¬ 
ing  the  modified  algorithm  on  image  5b  which  has  A/o  ■  1.  In  this  case 

using  the  above  analysis  procedures  we  chose  T  *  .83,  T  -  .6.  The  corre- 

11 

spondlng  error  probabilities  were  P£1  -  3.4xl0“4,  P£2  -  .76  and  P„,  -  .02. 
Figure  9b  is  the  post-filtered  version  of  9a. 


E3 
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VI.  EXTENSIONS 


Although  in  this  report  we  have  concentrated  on  the  simple  uniform 
Intensity  binary  segmentation  problem  we  point  out  that  the  approach  is 
easily  extended. 

First,  if  one  were  interested  in  the  colored  or  multi-spectral 

case,  one  could  model  the  image  as  a  set  of  intensity  matrices 

where  each  matrix  corresponded  to  the  intensity  level  of  a  primary  color 

or  of  a  given  spectral  band.  In  this  case  a  region  would  be  characterized 

by  the  set  of  intensities  {rf}^1 

k  1*1 

Next  observe  that  by  Increasing  the  number  of  allowable  Markov  states, 

the  multi-region  case  can  be  handled.  If  the  image  had  K  regions  then  a 

A  K 

segmentation  would  be  characterized  by  S  -  U  S,  . 

k-1  k 

Combining  these  two  extensions,  leads  to  the  following  generalization 


of  (5) 

A(S,{G£}) 

| S) 


MS)  +  !<{g£}|s) 


L  K 

C  -  E  [  E  (  E 
1-1  k-1  (i,j)eS 


-rk>  )] 


(34) 

(35) 


K 

*(S)  -EE  „  P...  (36) 

k-1  (i,j)eSk 

Using  (34) -(36)  estimation  algorithms  similar  to  those  described  earlier 
can  be  developed. 

Finally  we  point  out  that  this  approach  can  be  even  further  extended 
to  consider  textured  Images  where  region  textures  are  modelled  stochas¬ 
tically  using  probability  density  functions  or  dynamical  models  as  In 
[  8],  [ 9  J.  Zn  this  case  the  state  space  of  the  Markov  field  must  be 
modified.  That  is,  rather  than  characterize  a  region  by  a  simple  feature 

such  as  intensity,  one  could  define  it  to  be  a  parameter  vector  associated 
with  the  texture  model. 
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VII.  CONCLUDING  REMARKS 

This  paper  presented  a  new  segmentation  algorithm  based  upon  a 
simple  Markov  field  model  for  an  Image  and  MAP  estimation  techniques. 

The  two-dimensional  nature  of  the  image  forced  us  to  consider  a  subopt Imal 
algorithm  for  actual  segmentation.  We  also  presented  some  analytical 
techniques  for  estimating  key  algorithm  parameters.  Some  areas  for  future 
work  include  more  careful  study  of  the  extensions  outlined  in  the  previous 
section.  In  addition  it  would  be  of  Interest  to  explore  other  approxima¬ 
tions  to  the  true  MAP  estimate  and  compare  performance  of  these  alternative 
approaches  both  experimentally  and  analytically. 
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APPENDIX  A 


This  appendix  describes  the  dynamic  programming  techniques  necessary 
for  implementation  of  the  Stage  1  and  Stage  2  algorithms  discussed  in 
Section  III.  We  first  present  a  basic  forward  dynamic  programing  algo¬ 
rithm,  of  the  form  popularized  by  Bellman  [13],  for  segmenting  a  row  into 

J- 

regions,  §  ,  characterized  by  region  intensities  rR,  and  the  Gaussian 
noise  model.  Modification  of  this  algorithm  for  generating  the  M  most 
likely  segmentation  (the  Stage  1  algorithm)  are  then  presented.  Finally, 
the  Stage  2  algorithm  is  derived  by  simply  changing  a  few  of  the  defini¬ 
tions  introduced  for  row  segmentation. 

Forward  Dynamic  Programming  for  Segmentation  of  the  1th  Row  of  an  Image 
Preliminary  Definitions: 


'kj  ” 
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in(l-T) 

if  k*j 

N  &  N2 
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Set  c 
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max  +  Dfcc  +  it  )  ,  k»l,2,...K 

l<j<K  3  J 

maximizing  value  of  j,  k-l,...,K 


Step  3:  If  c<N  go  to  Step  2 
Step  4:  Set  £R(G(1,N2) ,S(i,N2)) 


max  l. 

!<)<*.  3 


p  -  maximizing  value  of  J 


Assign  pixel  <1,N2)  to  set  Sp(l,N2> 
Step  5:  c  •  c-1 


Assign  pixel  (i,c)  to  set  Sp(l,N2) 

Step  6:  If  c>l  go  to  5 
Step  7:  Stop. 

The  Stage  1  Algor it ha 

To  generate  the  M  most  likely  segmentations  modify  steps  2,  4,  and 
5  above  to 

Step  2':  m-max(i®“1,P  +  +  D.  ) ,  k-1 . K 

*  l£j«  3  K3  M-l . M 

l<p<M 
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Assign  pixel  (i,c)  to  set  S*  (i,N2)  ,  m-l,...,K  . 
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The  Stage  2  Algorithm 

The  Stage  2  algorithm  Is  a  simple  forward  dynamic  programming 
algorithm,  and  as  such,  can  be  obtained  from  the  basic  row  segmentation 
algorithm  by  redefining  the  following  parameters: 

Dkj  *  V6(J'N2)»*k(J*N2)) 

wkj  "  ic<§k(J,N2)) 


K  -  8 


In  addition  steps  4  and  5  should  be  redefined  as: 

Step  4":  t(G,S)  -  max  i*} 
l<j<K  3 

p •maximizing  value  of  j 
Assign  row  segmentation  S^CN.Nj)  to  § 
Step  5":  c  ■  c-1 


Assign  ^(c.N^)  to  S 


segmentation  with  T-.95  and  M-4. 


** 
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Figure  3a.  64x64  ellipse  in  additive  Gaussian  noise  Figure  3b.  Maximum  likelihood  threshold  segmentation, 

such  that  i/a  «  1. 


Figure  4a.  Intensity  profiles  for  correct  row  segmentation  S(1,N  )  and 
Incorrect  segmentation  S(i,N2>  for  case  1.  2 


s(i,n2) 


s(i,n2) 


Figure  4b.  Intensity  profiels  for  correct  row  segmentation  S(i,N^)  and 
incorrect  segmentation  S(i,N2)  for  case  2.  2 


Figure  5e.  MAP  segmentation  with  T*.3  and  K*8.  Figure  5f .  MAP  segmentation  with  1 


Figure  6c.  MAP  segmentation  with  T*.95  and  M»8. 
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